function [MRPSO,strategy] = MRPSO(task,execution_time,task_num,task_out,task_in,LOPSO,MaxNum)
%% 初始值设置(粒子群算法参数）
particlesize = 100;                     % 粒子群规模
c1 = 2;                                 % 个体经验学习因子
c2 = 2;
c3 = 2;
c4 = 2;
c5 = 2;
c6 = 2;
alpha = 0.9;                            % 降温系数，解决局部最优
wmax=0.9;
wmin=0.4;
narvs=task_num+7;

%% 粒子群算法——初始位置、初始fitness值计算
f=zeros(1,particlesize);                % 预留空间
vmax = 0.5;                             % 粒子的最大飞行速度
vvmax= 0.5;                             % 订购方案最大飞行速度
x =  rand(particlesize, narvs);         % 粒子所在的位置
v =0.2*(rand(particlesize,narvs)-0.5);  % 粒子的飞翔速度  生成每个粒子的飞翔速度，由于是只有一个变量，所以速度是一维的
temp=to_ceil(x,task_num);
for i= 1:particlesize
    f(i) = object_value(temp(i,:),task,execution_time,task_num,task_out,task_in);
end
[sort_f,i]=sort(f);
p=[sort_f(1:3);i(1:3)];                     %找到最好的三个
best=x(p(2,1),:);
best_value=p(1,1);                          %群体最优
best_array(1)= best_value;
f_best=f;                                   %个体最优适应度值
f_best_pop=x;                               %个体最优位置
T=LOPSO;                                    %设置初始温度
Best_value_all=best_value;

%% 粒子群算法——迭代寻优
xx=zeros(particlesize,narvs);
time=1;
time_num=MaxNum;
while(time~=time_num)                         % 开始迭代
    w =wmax-((wmax-wmin)/time_num)*time;                                 % 动态惯性因子
    for k=1:particlesize
        %更新粒子速度
        %v(k,:) = w*v(k,:) + c1*rand*(f_best_pop(k,:) - x(k,:)) + c2*rand*(best - x(k,:));
        %v(k,1:7)=w*v(k,1:7)+c1*rand*(x(p(2,1),1:7)-xx(k,1:7))+c2*rand*(x(p(2,2),1:7)-xx(k,1:7))+c3*rand*(x(p(2,3),1:7)-xx(k,1:7))+c4*rand*(x(k,1:7)-xx(k,1:7))+c5*(rand(1,7)-0.5);
        %v(k,8:narvs)=w*v(k,8:narvs)+c1*rand*(x(p(2,1),8:narvs)-xx(k,8:narvs))+c2*rand*(x(p(2,2),8:narvs)-xx(k,8:narvs))+c3*rand*(x(p(2,3),8:narvs)-xx(k,8:narvs))+c4*rand*(x(k,8:narvs)-xx(k,8:narvs))+c6*(rand(1,narvs-7)-0.5);
        v(k,1:7)=w*v(k,1:7)+0.1*(c1*rand*(x(p(2,1),1:7)-xx(k,1:7))+c2*rand*(x(p(2,2),1:7)-xx(k,1:7))+c3*rand*(x(p(2,3),1:7)-xx(k,1:7))+c4*rand*(x(k,1:7)-xx(k,1:7))+c5*(rand(1,7)-0.5));
        v(k,8:narvs)=w*v(k,8:narvs)+0.1*(c1*rand*(x(p(2,1),8:narvs)-xx(k,8:narvs))+c2*rand*(x(p(2,2),8:narvs)-xx(k,8:narvs))+c3*rand*(x(p(2,3),8:narvs)-xx(k,8:narvs))+c4*rand*(x(k,8:narvs)-xx(k,8:narvs))+c6*(rand(1,narvs-7)-0.5));
        %粒子速度越界处理，粒子位置越界处理
        for i=1 : 7
            if v(k,i)>vvmax
                v(k,i)=vvmax;
            end
            if v(k,i)<-vvmax
                v(k,i)=-vvmax;
            end
            %更新粒子位置
            xx(k,i)=x(k,i)+v(k,i);
            if xx(k,i)>1
                xx(k,i)=rand;
            elseif xx(k,i)<0
                xx(k,i)=rand;
            end
        end
        for i=8 : narvs
            if v(k,i)>vmax
                v(k,i)=vmax;
            end
            if v(k,i)<-vmax
                v(k,i)=-vmax;
            end
            %更新粒子位置
            xx(k,i)=x(k,i)+v(k,i);
            if xx(k,i)>1
                xx(k,i)=1-10e-10;
            elseif xx(k,i)<0
                xx(k,i)=10e-10;
            end
        end
    end
    temp=to_ceil(xx,task_num);              %x处理成串
    %更新适应度值
    for i=1:particlesize
        f(i) = object_value(temp(i,:),task,execution_time,task_num,task_out,task_in);
    end
    for i=1:particlesize                    %找到局部最优
        %更新个体最优
        if f(i)<f_best(i)                   %如果这个粒子此次迭代更好
            f_best(i)=f(i);                 %适应度值
            f_best_pop(i,:)=xx(i,:);        %位置
            x(i,:)=xx(i,:);                 %更新个体最优位置
        elseif(rand < exp(-(f(i)-f_best(i))/T)) %以一定概率接收差的迭代结果
            f_best(i)=f(i);                 %适应度值
            f_best_pop(i,:)=xx(i,:);        %位置
            x(i,:)=xx(i,:);
        end
    end
    %更新群体最优
    [sort_f,i]=sort(f_best);
    pp=[sort_f(1:3);i(1:3)];                 %找到最好的三个
    best1=x(pp(2,1),:);                      %此次迭代最好的粒子
    best_value1=pp(1,1);
    if(best_value1<best_value)
        best=best1;                         %位置
        best_value=best_value1;             %适应度值
        p=pp;
    elseif(rand < exp(-(best_value1-best_value)/T))
        best=best1;                         %位置
        best_value=best_value1;             %适应度值
        p=pp;
    end
    
    if(Best_value_all>best_value)           %历史最佳值
        Best_value_all=best_value;
    end
    best_array(time+1) =best_value;
    T = alpha*T;                            %降温
    time=time+1;
    if(time==time_num)   %最后一次迭代
        if(best_value>Best_value_all)       %若迭代结果最优差于历史最优，重新进行迭代
            disp("迭代结果最优差于历史最优，重新进行迭代");
            time_num=time_num+MaxNum;
        end
        T=LOPSO/alpha;            %设置新的初始温度
        alpha=alpha*0.9;
        if(alpha<0.5)
            disp("迭代次数过多且尚未得到最优迭代值，重新迭代");
            %% 粒子群算法——初始位置、初始fitness值计算
            f=zeros(1,particlesize);                % 预留空间
            w =0.8;                                 % 惯性因子
            vmax = 0.05;                            % 粒子的最大飞翔速度
            x =  rand(particlesize, narvs);         % 粒子所在的位置
            v =0.2*(rand(particlesize,narvs)-0.5);  % 粒子的飞翔速度  生成每个粒子的飞翔速度，由于是只有一个变量，所以速度是一维的
            temp=to_ceil(x,task_num);
            for i= 1:particlesize
                f(i) = object_value(temp(i,:),task,execution_time,task_num,task_out,task_in);
            end
            [sort_f,i]=sort(f);
            p=[sort_f(1:3);i(1:3)];                     %找到最好的三个
            best=x(p(2,1),:);
            best_value=p(1,1);                          %群体最优
            f_best=f;                                   %个体最优适应度值
            f_best_pop=x;                               %个体最优位置
            T=LOPSO;                                       %设置初始温度
            if(Best_value_all>best_value)               %历史最佳值
                Best_value_all=best_value;
            end
            alpha=0.9;
        end
    end
end

%% 粒子群算法——输出结果
best=to_ceil(best,task_num);
if best(1)<3                                %best修正
    for b=best(1)+2:4
        best(b)=0;
    end
    for c=best(1)+5:7
        best(c)=0;
    end
end
disp("最佳配置");best
strategy=best;
figure();
plot(best_array);
title("MRPSO");
MRPSO=best_value;
end

